#include "AnaDef.h"
#include "XGLUtils.h"
#include "AliTRDdEdxBaseUtils.h"

void doSeparation(const Int_t itype, const Bool_t kinvq, const Int_t ib, const Int_t ic)
{
  printf("doSeparation: itype %d, kinvq %d, ib %d, ic %d\n", itype, kinvq, ib, ic);

  gErrorIgnoreLevel = 1001;
  
  TTree * tree = (TTree*)gDirectory->Get("newdraw_Tree");
  TString hname("hh");
  TString cut("1");
  TString var;
  TH2D *hh=0x0;

  if(kinvq){  
    tree->SetAlias("nch", "trdnch1");
    tree->SetAlias("ncls", "trdncls1");
    tree->SetAlias("trdq", "trdQ1");
  }
  else{
    tree->SetAlias("nch", "trdnch0");
    tree->SetAlias("ncls", "trdncls0");
    tree->SetAlias("trdq", "trdQ0");
  }

  tree->SetAlias("rtpc","tpcsig/AliTRDdEdxBaseUtils::QMeanTPC(pow(10,ltbg))");

  //==============================================

  hname+=Form("type%dQ%db%dc%d",itype, kinvq,ib,ic);

  //==============================================

  if(itype>=0 && itype<=1){//cosmic
    cut+= Form("&& (abs(mag)==1 || itype==%d)",itype);

    cut+= "&& !(abs(mag)==1 && ltbg>1.4)";
    cut+= "&& !(abs(mag)==5 && ltbg<1.5)";

    tree->SetAlias("ele1kg","abs(mag)==1&&tpcsig>44&&ltbg<1.5");
    tree->SetAlias("elefeb","abs(mag)==5 && run<150000 && ((tpcsig>160&&ltbg<1.5)||(tpcsig>59&&tpcsig<90&&ltbg<1.6))");
    cut+= "&& !elefeb && !ele1kg";
  }
  else if(itype==2){//pp
    cut+= Form("&& itype==%d",itype);

    //cut+= Form("&& rtpc>0.8 && rtpc<1.2");
    //cut+= Form("&& ncls>115");
    printf("no ncls cut yet!\n");

    //printf("************** 10 d + e !!!!!!!\n");
    //cut+= "&& run >= 121527 && run <= 126460"; printf("selecting 10ddddddddddd !!\n");
    //cut+= "&& run >= 126461  && run <= 130930"; printf("selecting 10eeeeeeeeee !!\n");
  }
  else if(itype==3){//Pb all centrality
    cut+= Form("&& itype>=10");

    cut+= Form("&& rtpc>0.95 && rtpc<1.3");
  }
  else if(itype==4){//Pb most central
    cut+= Form("&& itype>=10 && itype<=11");
  }
  else{
    printf("wrong itype! %d\n", itype); exit(1);
  }

  if(ib==0)      cut+= "&& mag <0";
  else if(ib==1) cut+= "&& mag >0";

  if(ic==0)      cut+= "&& charge <0";
  else if(ic==1) cut+= "&& charge >0";

  cut+= "&& nch==6"; 

  //========================================================

  var = "trdq:ltbg"; hh=new TH2D(hname,"", 50, -0.5, 4.5, 20000, 0, 5);

  // var = "tpcsig:ltbg"; hh=new TH2D(hname+"tpc","", 50, -0.5, 4.5, 10000, 0, 200);
  
  //==================================

  cout<<hname<<endl;
  cout<<var<<endl;                    
  cout<<cut<<endl;
  tree->SetScanField(0);
  tree->Scan(var,cut,"",50);
  
  const Int_t totn = tree->Draw(Form("%s>>%s", var.Data(), hh->GetName()),cut);
  if(totn<100){
    printf("too few events! %d\n", totn);
    return;
  }

  TH1D *hfit = new TH1D("hfit","", 500, -0.5, 4.5);
  for(Int_t ii=1; ii<=hfit->GetNbinsX(); ii++){
    Double_t xx = hfit->GetBinCenter(ii);
    Double_t yy = AliTRDdEdxBaseUtils::Q0MeanTRDpp(pow(10,xx));
    hfit->SetBinContent(ii,yy);
  }
  //hfit->SetLineColor(kRed);
  hfit->Draw("same");
  //=====================================
  
  TH1D *hnor=0, *hmpv=0, *hwid=0, *hres=0, *hchi=0;

  XGLUtils::FitSlicesY(hh, 0x0, hnor, hmpv, hwid, hres, hchi, 50, 0, 0.99);//if highFrac=1, then hwid fluctuated a lot!

  const Int_t anchor = hmpv->FindBin(0.55);
  const Int_t target = hmpv->FindBin(3.15);
  const Double_t xanc = hmpv->GetBinCenter(anchor);
  const Double_t xtar = hmpv->GetBinCenter(target);

  const Double_t scale = hmpv->GetBinContent(anchor);

  //=====================================
 
  if( var.Contains("trdq") ){
     hmpv->Scale(1/scale);
    //very important also scale hwid
    hwid->Scale(1/scale);
  }

  //=====================================

  const Double_t a0 = hnor->GetBinContent(anchor);
  const Double_t a1 = hnor->GetBinContent(target);
  const Double_t m0 = hmpv->GetBinContent(anchor);
  const Double_t m1 = hmpv->GetBinContent(target);
  const Double_t w0 = hwid->GetBinContent(anchor);
  const Double_t w1 = hwid->GetBinContent(target);
  const Double_t guess = fabs(m1-m0)/(w0+w1)*2;
  
  printf("%s TotN %d , Scale@ %.2f : %.4f , Sep@ %.2f : %.4f , N %.0f %.0f , Mean %.1f %.4f , Wid %.4f  %.4f , Res %.4f %.4f\n", hh->GetName(), totn, xanc, scale, xtar, guess, a0, a1, m0, m1, w0, w1, w0/m0, w1/m1);

  TH1D *hsep=0x0;
  hsep= XGLUtils::GetSeparation(hmpv, hwid, anchor, kFALSE);

  //=================================
  if(hh->GetNbinsY()>1000) hh->RebinY(100);
  TH2D *h2=XGLUtils::NormalHist(hh, 50,1);

  TFile *fout=new TFile(Form("out%s.root",hname.Data()),"recreate");
  hh->Write();
  h2->Write();
  hnor->Write();
  hmpv->Write();
  hwid->Write();
  hres->Write();
  hchi->Write();
  hsep->Write();
  hfit->Write();
  fout->Save();
  fout->Close();  
}

int main(const int argc, char *argv[])
{
  Int_t nums[argc];
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
    nums[ii] = atoi(argv[ii]);
  }

  if(argc!=6){
    printf("argc!=5\n");
    return 1;
  }

  TFile::Open(argv[1]);
  doSeparation(nums[2], nums[3], nums[4], nums[5]);
  return 0;
}
